This file is used to analyse the immune cells dataset.
library(dplyr)
library(patchwork)
library(ggplot2)
library(ComplexHeatmap)
library(org.Mm.eg.db)
.libPaths()
## [1] "/usr/local/lib/R/library"
In this section, we set the global settings of the analysis. We will store data there :
save_name = "iblmors"
out_dir = "."
We load the dataset :
sobj = readRDS(paste0(out_dir, "/", save_name, "_sobj.rds"))
sobj
## An object of class Seurat
## 16701 features across 3532 samples within 1 assay
## Active assay: RNA (16701 features, 2000 variable features)
## 6 dimensional reductions calculated: RNA_pca, RNA_pca_20_tsne, RNA_pca_20_umap, harmony, harmony_20_umap, harmony_20_tsne
We load the sample information :
sample_info = readRDS(paste0(out_dir, "/../../1_metadata/hs_hd_sample_info.rds"))
project_names_oi = sample_info$project_name
graphics::pie(rep(1, nrow(sample_info)),
col = sample_info$color,
labels = sample_info$project_name)
Here are custom colors for each cell type :
color_markers = readRDS(paste0(out_dir, "/../../1_metadata/hs_hd_color_markers.rds"))
data.frame(cell_type = names(color_markers),
color = unlist(color_markers)) %>%
ggplot2::ggplot(., aes(x = cell_type, y = 0, fill = cell_type)) +
ggplot2::geom_point(pch = 21, size = 5) +
ggplot2::scale_fill_manual(values = unlist(color_markers), breaks = names(color_markers)) +
ggplot2::theme_classic() +
ggplot2::theme(legend.position = "none",
axis.line = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank(),
axis.text.y = element_blank(),
axis.text.x = element_text(angle = 30, hjust = 1))
This is the projection of interest :
name2D = "harmony_20_tsne"
We design a custom function to make the GSEA plot and a word cloud graph :
make_gsea_plot = function(gsea_results, gs_oi, fold_change, metric = "FC") {
fold_change$metric = fold_change[, metric]
plot_list = lapply(gs_oi, FUN = function(gene_set) {
# Gene set content
gs_content = gene_sets %>%
dplyr::filter(gs_name == gene_set) %>%
dplyr::pull(ensembl_gene) %>%
unique()
# Gene set size
nb_genes = length(gs_content)
# Enrichment metrics
NES = gsea_results@result[gene_set, "NES"]
p.adjust = gsea_results@result[gene_set, "p.adjust"] %>%
round(., 4)
qvalues = gsea_results@result[gene_set, "qvalues"]
if (p.adjust > 0.05) {
p.adjust = paste0("<span style='color:red;'>", p.adjust, "</span>")
}
my_subtitle = paste0("\nNES : ", round(NES, 2),
" | padj : ", p.adjust,
" | qval : ", round(qvalues, 4),
" | set size : ", nb_genes, " genes")
# Size limits
lower_FC = min(fold_change[gs_content, ]$metric, na.rm = TRUE)
upper_FC = max(fold_change[gs_content, ]$metric, na.rm = TRUE)
# Plot
p = enrichplot::gseaplot2(x = gsea_results, geneSetID = gene_set) +
ggplot2::labs(title = gene_set,
subtitle = my_subtitle) +
ggplot2::theme(plot.title = element_text(hjust = 0.5, face = "bold",
margin = ggplot2::margin(3, 3, 5, 3)),
plot.subtitle = ggtext::element_markdown(hjust = 0.5,
size = 10))
wc = ggplot2::ggplot(fold_change[gs_content, ],
aes(label = gene_name, size = abs(metric), color = metric)) +
ggwordcloud::geom_text_wordcloud_area(show.legend = TRUE) +
ggplot2::scale_color_gradient2(
name = metric,
low = aquarius::color_cnv[1],
mid = "gray70", midpoint = 0,
high = aquarius::color_cnv[3]) +
ggplot2::scale_size_area(max_size = 7) +
ggplot2::theme_minimal() +
ggplot2::guides(size = "none")
return(list(p, wc))
}) %>% unlist(., recursive = FALSE)
return(plot_list)
}
We visualize gene expression for some markers :
features = c("percent.mt", "percent.rb", "nFeature_RNA")
plot_list = lapply(features, FUN = function(one_gene) {
Seurat::FeaturePlot(sobj, features = one_gene,
reduction = name2D) +
ggplot2::theme(aspect.ratio = 1) +
ggplot2::scale_color_gradientn(colors = aquarius::color_gene) +
Seurat::NoAxes()
})
patchwork::wrap_plots(plot_list, ncol = 3)
We visualize clusters :
cluster_plot = Seurat::DimPlot(sobj, reduction = name2D, label = TRUE) +
Seurat::NoAxes() +
ggplot2::theme(aspect.ratio = 1)
cluster_plot
We visualize cell type split by sample :
plot_list = aquarius::plot_split_dimred(sobj,
reduction = name2D,
split_by = "project_name",
group_by = "cell_type",
split_color = setNames(sample_info$color,
nm = sample_info$project_name),
group_color = color_markers,
bg_pt_size = 0.5, main_pt_size = 0.5)
plot_list[[length(plot_list) + 1]] = cluster_plot
patchwork::wrap_plots(plot_list, ncol = 4) &
Seurat::NoLegend()
We summarize major cell type by cluster :
cell_type_clusters = sobj@meta.data[, c("cell_type", "seurat_clusters")] %>%
table() %>%
prop.table(., margin = 2) %>%
apply(., 2, which.max)
cell_type_clusters = setNames(levels(sobj$cell_type)[cell_type_clusters],
nm = names(cell_type_clusters))
We define cluster type :
sobj$cluster_type = cell_type_clusters[sobj$seurat_clusters] %>%
as.factor()
table(sobj$cluster_type, sobj$cell_type)
##
## CD4 T cells CD8 T cells Langerhans cells macrophages B cells cuticle
## IBL 0 0 0 0 1 1
## ORS 2 0 0 0 3 7
##
## cortex medulla IRS proliferative IBL ORS IFE HFSC sebocytes
## IBL 0 0 0 9 1465 5 7 19 3
## ORS 2 20 7 37 26 1785 65 48 20
We subset color_markers :
color_markers = color_markers[levels(sobj$cluster_type)]
We compare cluster annotation and cell type annotation :
p1 = Seurat::DimPlot(sobj, group.by = "cell_type",
reduction = name2D, cols = color_markers) +
ggplot2::labs(title = "Cell type") +
Seurat::NoAxes() +
ggplot2::theme(aspect.ratio = 1,
plot.title = element_text(hjust = 0.5))
p2 = Seurat::DimPlot(sobj, group.by = "cluster_type",
reduction = name2D, cols = color_markers) +
ggplot2::labs(title = "Cluster type") +
Seurat::NoAxes() +
ggplot2::theme(aspect.ratio = 1,
plot.title = element_text(hjust = 0.5))
patchwork::wrap_plots(p1, p2, guides = "collect")
We make a barplot to compare IBL and ORS populations in HS vs HD samples. The proportion of IBL is indicated on top of bars.
quantif = table(sobj$sample_identifier,
sobj$cluster_type) %>%
as.data.frame.table() %>%
`colnames<-`(c("Sample", "cell_type", "nb_cells")) %>%
dplyr::group_by(Sample) %>%
dplyr::mutate(total_cells = sum(nb_cells)) %>%
as.data.frame() %>%
dplyr::filter(cell_type == "IBL") %>%
dplyr::mutate(prop_ibl = nb_cells / total_cells) %>%
dplyr::mutate(prop_ibl = 100*round(prop_ibl, 4))
sobj$seurat_clusters = factor(sobj$seurat_clusters,
levels = names(sort(cell_type_clusters)))
aquarius::plot_barplot(df = table(sobj$sample_identifier,
sobj$seurat_clusters) %>%
as.data.frame.table() %>%
`colnames<-`(c("sample_identifier", "clusters", "nb_cells")),
x = "sample_identifier", y = "nb_cells", fill = "clusters",
position = position_fill()) +
ggplot2::scale_fill_manual(values = c(colorRampPalette(c("chartreuse1", "chartreuse4"))(table(sort(cell_type_clusters))["IBL"]),
colorRampPalette(c("cadetblue1", "cadetblue4"))(table(sort(cell_type_clusters))["ORS"])),
breaks = names(sort(cell_type_clusters)),
name = "Cell type") +
ggplot2::geom_label(data = quantif, inherit.aes = FALSE,
aes(x = .data$Sample, y = 1.05, label = .data$prop_ibl),
label.size = 0, size = 5)
In this section, we perform DE between :
We save the results in a list :
list_results = list()
We make over-representation analysis for each group of genes. We load gene sets from MSigDB :
gene_sets = aquarius::get_gene_sets(species = "Homo sapiens")
gene_sets = gene_sets$gene_sets
head(gene_sets)
## # A tibble: 6 x 16
## gs_cat gs_subcat gs_name gene_symbol entrez_gene ensembl_gene human_gene_symb~
## <chr> <chr> <chr> <chr> <int> <chr> <chr>
## 1 C5 GO:BP GOBP_1~ AASDHPPT 60496 ENSG0000014~ AASDHPPT
## 2 C5 GO:BP GOBP_1~ ALDH1L1 10840 ENSG0000014~ ALDH1L1
## 3 C5 GO:BP GOBP_1~ ALDH1L2 160428 ENSG0000013~ ALDH1L2
## 4 C5 GO:BP GOBP_1~ MTHFD1 4522 ENSG0000010~ MTHFD1
## 5 C5 GO:BP GOBP_1~ MTHFD1L 25902 ENSG0000012~ MTHFD1L
## 6 C5 GO:BP GOBP_1~ MTHFD2L 441024 ENSG0000016~ MTHFD2L
## # ... with 9 more variables: human_entrez_gene <int>, human_ensembl_gene <chr>,
## # gs_id <chr>, gs_pmid <chr>, gs_geoid <chr>, gs_exact_source <chr>,
## # gs_url <chr>, gs_description <chr>, category <chr>
How many gene sets ?
gene_sets[, c("gs_subcat", "gs_name")] %>%
dplyr::distinct() %>%
dplyr::pull(gs_subcat) %>%
table() %>%
as.data.frame.table() %>%
`colnames<-`(c("Category", "Nb gene sets"))
## Category Nb gene sets
## 1 50
## 2 CP:KEGG 186
## 3 CP:PID 196
## 4 CP:REACTOME 1615
## 5 CP:WIKIPATHWAYS 664
## 6 GO:BP 7658
## 7 GO:CC 1006
## 8 GO:MF 1738
We get gene name and gene ID correspondence :
gene_corresp = sobj@assays[["RNA"]]@meta.features[, c("gene_name", "Ensembl_ID")] %>%
`colnames<-`(c("NAME", "ID")) %>%
dplyr::mutate(ID = as.character(ID))
rownames(gene_corresp) = gene_corresp$ID
head(gene_corresp)
## NAME ID
## ENSG00000238009 AL627309.1 ENSG00000238009
## ENSG00000237491 AL669831.5 ENSG00000237491
## ENSG00000225880 LINC00115 ENSG00000225880
## ENSG00000230368 FAM41C ENSG00000230368
## ENSG00000230699 AL645608.3 ENSG00000230699
## ENSG00000187634 SAMD11 ENSG00000187634
group_name = "IBL_vs_ORS"
We change cell identities to cluster type :
Seurat::Idents(sobj) = sobj$cluster_type
table(Seurat::Idents(sobj), sobj$sample_type)
##
## HS HD
## IBL 1284 226
## ORS 1563 459
We identify specific markers for each population :
mark = Seurat::FindMarkers(sobj, ident.1 = "IBL", ident.2 = "ORS")
mark = mark %>%
dplyr::filter(p_val_adj < 0.05) %>%
dplyr::arrange(-avg_logFC, pct.1 - pct.2)
list_results[[group_name]]$mark = mark
dim(mark)
## [1] 1027 5
head(mark, n = 20)
## p_val avg_logFC pct.1 pct.2 p_val_adj
## KRT16 0.000000e+00 4.323925 0.920 0.116 0.000000e+00
## KRT6B 0.000000e+00 3.627468 0.921 0.285 0.000000e+00
## FABP5 0.000000e+00 3.348761 0.993 0.509 0.000000e+00
## KRT17 0.000000e+00 3.020736 0.987 0.733 0.000000e+00
## KRT6A 1.267920e-155 2.994730 0.764 0.476 2.117554e-151
## CST6 1.834243e-160 2.741295 0.391 0.038 3.063369e-156
## KRT6C 9.379153e-283 2.653075 0.568 0.038 1.566412e-278
## TMSB4X 0.000000e+00 2.356496 0.989 0.864 0.000000e+00
## LGALS1 0.000000e+00 2.320083 0.929 0.484 0.000000e+00
## S100A2 0.000000e+00 2.134213 0.996 0.977 0.000000e+00
## CALML3 0.000000e+00 2.003586 0.972 0.686 0.000000e+00
## GJA1 0.000000e+00 1.973282 0.875 0.586 0.000000e+00
## GJB6 0.000000e+00 1.907928 0.916 0.656 0.000000e+00
## SBSN 9.451796e-242 1.751399 0.534 0.047 1.578544e-237
## APOE 0.000000e+00 1.744652 0.983 0.772 0.000000e+00
## NDUFA4L2 1.058296e-239 1.730465 0.697 0.217 1.767461e-235
## LYPD3 0.000000e+00 1.680352 0.738 0.135 0.000000e+00
## TUBB2A 0.000000e+00 1.667800 0.756 0.181 0.000000e+00
## SDC1 0.000000e+00 1.662321 0.838 0.517 0.000000e+00
## GJB2 1.288894e-254 1.634612 0.826 0.595 2.152582e-250
There are 1027 genes differentially expressed. How many for each population ?
# IBL
mark_ibl = mark %>%
dplyr::filter(avg_logFC > 0)
nrow(mark_ibl)
## [1] 434
# ORS
mark_ors = mark %>%
dplyr::filter(avg_logFC < 0)
nrow(mark_ors)
## [1] 593
We represent the information on a figure :
mark$gene_name = rownames(mark)
mark_to_label = rbind(
# up-regulated in IBL
mark %>% dplyr::top_n(., n = 20, wt = avg_logFC),
# up-regulated in ORS
mark %>% dplyr::top_n(., n = 20, wt = -avg_logFC),
# representative and selective for IBL
mark %>% dplyr::top_n(., n = 20, wt = (pct.1 - pct.2)),
# representative and selective for ORS
mark %>% dplyr::top_n(., n = 20, wt = -(pct.1 - pct.2))) %>%
dplyr::distinct()
ggplot2::ggplot(mark, aes(x = pct.1, y = pct.2, col = avg_logFC)) +
ggplot2::geom_abline(slope = 1, intercept = 0, lty = 2) +
ggplot2::geom_point() +
ggrepel::geom_label_repel(data = mark_to_label, max.overlaps = Inf,
aes(x = pct.1, y = pct.2, label = gene_name),
col = "black", fill = NA, size = 3, label.size = NA) +
ggplot2::labs(title = "Differentially expressed genes",
subtitle = "between IBL and ORS") +
ggplot2::scale_color_gradient2(low = aquarius::color_cnv[1],
mid = aquarius::color_cnv[2],
high = aquarius::color_cnv[3],
midpoint = 0) +
ggplot2::theme_classic() +
ggplot2::theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
We represent the best genes on a violin plot, group by cluster type and split by sample type :
mark_to_label = mark_to_label %>%
dplyr::arrange(pct.1 - pct.2)
plot_list = lapply(mark_to_label$gene_name, FUN = function(gene) {
p = Seurat::VlnPlot(sobj, features = gene, pt.size = 0.001,
group.by = "cluster_type",
split.by = "sample_type") +
ggplot2::scale_fill_manual(breaks = c("HS", "HD"),
values = c("#C55F40", "#2C78E6")) +
ggplot2::theme(axis.title.x = element_blank(),
axis.title.y = element_blank(),
legend.position = "none")
return(p)
})
patchwork::wrap_plots(plot_list, ncol = 5)
On the figure, we can see that some specific markers for a population are indeed specific to a sample type rather than the whole population.
We explore enrichment in gene sets for IBL population
genes_of_interest = rownames(mark_ibl)
enrichr_results = aquarius::run_enrichr(gene_names = genes_of_interest,
gene_corresp = gene_corresp,
gene_sets = gene_sets[, c("gs_name", "ensembl_gene")],
make_plot = TRUE,
plot_title = "Up-regulated in IBL compared to ORS")
list_results[[group_name]]$enrichr_ibl = enrichr_results$ego
enrichr_results$plot +
ggplot2::theme(axis.text.y = element_text(size = 8))
We explore enrichment in gene sets for ORS population.
genes_of_interest = rownames(mark_ors)
enrichr_results = aquarius::run_enrichr(gene_names = genes_of_interest,
gene_corresp = gene_corresp,
gene_sets = gene_sets[, c("gs_name", "ensembl_gene")],
make_plot = TRUE,
plot_title = "Up-regulated in ORS compared to IBL")
list_results[[group_name]]$enrichr_ors = enrichr_results$ego
enrichr_results$plot +
ggplot2::theme(axis.text.y = element_text(size = 8))
We run a GSEA for all gene sets, from the full count matrix :
ranked_gene_list = aquarius::run_foldchange(Seurat::GetAssayData(sobj, assay = "RNA", slot = "counts"),
group1 = colnames(sobj)[sobj@active.ident %in% "IBL"],
group2 = colnames(sobj)[sobj@active.ident %in% "ORS"])
names(ranked_gene_list) = gene_corresp$ID
gsea_results = aquarius::gsea_run(ranked_gene_list = ranked_gene_list,
gene_sets = gene_sets[, c("gs_name", "ensembl_gene")],
GSEA_p_val_thresh = 1)
list_results[[group_name]]$gsea = gsea_results
gsea_results@result %>%
dplyr::filter(pvalue < 0.05) %>%
dplyr::top_n(., n = 200, wt = abs(NES)) %>%
dplyr::mutate(too_long = ifelse(nchar(ID) > 60, yes = TRUE, no = FALSE)) %>%
dplyr::mutate(ID = stringr::str_sub(ID, end = 60)) %>%
dplyr::mutate(ID = ifelse(too_long, yes = paste0(ID, "..."), no = ID)) %>%
aquarius::gsea_plot(show_legend = TRUE) +
ggplot2::labs(title = "GSEA using all genes (count matrix)") +
ggplot2::theme(plot.title = element_text(size = 20))
We make the gsea plot for two gene sets :
p1 = enrichplot::gseaplot2(x = gsea_results, geneSetID = "REACTOME_KERATINIZATION") +
ggplot2::labs(title = "REACTOME_KERATINIZATION") +
ggplot2::theme(plot.title = element_text(hjust = 0.5, face = "bold",
margin = ggplot2::margin(3, 3, 5, 3)))
p2 = enrichplot::gseaplot2(x = gsea_results, geneSetID = "HALLMARK_INTERFERON_GAMMA_RESPONSE") +
ggplot2::labs(title = "HALLMARK_INTERFERON_GAMMA_RESPONSE") +
ggplot2::theme(plot.title = element_text(hjust = 0.5, face = "bold",
margin = ggplot2::margin(3, 3, 5, 3)))
p1 | p2
Maybe those gene sets are specific to sample type rather than the whole cell population.
group_name = "cluster5_vs_ORS"
We subset the Seurat object and change cell identities to sample type.
subsobj = subset(sobj, cluster_type == "ORS")
subsobj$seurat_clusters = base::droplevels(subsobj$seurat_clusters)
Seurat::Idents(subsobj) = subsobj$seurat_clusters
table(subsobj$seurat_clusters)
##
## 0 1 3 5 7
## 669 438 384 302 229
We identify specific markers for each population :
mark = Seurat::FindMarkers(subsobj, ident.1 = 5)
mark = mark %>%
dplyr::filter(p_val_adj < 0.05) %>%
dplyr::arrange(-avg_logFC, pct.1 - pct.2)
list_results[[group_name]]$mark = mark
dim(mark)
## [1] 395 5
head(mark, n = 20)
## p_val avg_logFC pct.1 pct.2 p_val_adj
## SLPI 7.619364e-128 1.5306935 0.765 0.164 1.272510e-123
## MT1X 4.989262e-73 1.4231367 0.954 0.680 8.332567e-69
## KRT15 2.275685e-123 1.3705217 1.000 0.874 3.800622e-119
## LY6D 1.984327e-88 1.2250993 0.503 0.080 3.314025e-84
## IGFBP3 1.494292e-140 1.1855371 0.437 0.016 2.495618e-136
## IFI27 1.438175e-99 1.1442636 0.556 0.085 2.401895e-95
## MT1E 2.040603e-70 1.0839955 0.957 0.659 3.408012e-66
## WNT3 1.095829e-212 1.0819468 0.924 0.134 1.830143e-208
## AQP3 5.658795e-83 1.0173107 0.977 0.733 9.450753e-79
## NEAT1 3.110435e-52 0.9735049 0.990 0.905 5.194737e-48
## TSC22D3 1.298633e-75 0.9601088 0.742 0.258 2.168848e-71
## MTRNR2L1 4.400300e-32 0.9329154 0.815 0.593 7.348942e-28
## CXCL14 7.450833e-63 0.9282610 0.997 0.840 1.244364e-58
## ZFP36 1.579077e-42 0.9002595 0.742 0.406 2.637216e-38
## GLUL 7.969325e-58 0.8930699 0.934 0.699 1.330957e-53
## ID1 2.365404e-17 0.8903797 0.652 0.442 3.950461e-13
## FOS 1.096048e-36 0.8789679 0.967 0.837 1.830511e-32
## IL1R2 5.014612e-122 0.8759829 0.675 0.110 8.374904e-118
## IL18 6.622739e-101 0.8511305 0.970 0.456 1.106064e-96
## HOPX 1.939070e-75 0.8395659 0.871 0.417 3.238441e-71
There are 395 genes differentially expressed. How many for each population ?
# Cluster 5
mark_cluster5 = mark %>%
dplyr::filter(avg_logFC > 0)
nrow(mark_cluster5)
## [1] 225
# Other ORS
mark_not5 = mark %>%
dplyr::filter(avg_logFC < 0)
nrow(mark_not5)
## [1] 170
We represent the information on a figure :
mark$gene_name = rownames(mark)
mark_signif = rbind(
# up-regulated in cluster 5
mark %>% dplyr::top_n(., n = 20, wt = avg_logFC),
# up-regulated in other ORS
mark %>% dplyr::top_n(., n = 20, wt = -avg_logFC),
# representative and selective for cluster 5
mark %>% dplyr::top_n(., n = 20, wt = (pct.1 - pct.2)),
# representative and selective for other ORS
mark %>% dplyr::top_n(., n = 20, wt = -(pct.1 - pct.2))) %>%
dplyr::distinct()
ggplot2::ggplot(mark, aes(x = pct.1, y = pct.2, col = avg_logFC)) +
ggplot2::geom_abline(slope = 1, intercept = 0, lty = 2) +
ggplot2::geom_point() +
ggrepel::geom_label_repel(data = mark_signif, max.overlaps = Inf,
aes(x = pct.1, y = pct.2, label = gene_name),
col = "black", fill = NA, size = 3, label.size = NA) +
ggplot2::labs(title = "Differentially expressed genes",
subtitle = "between cluster 5 and other ORS") +
ggplot2::scale_color_gradient2(low = aquarius::color_cnv[1],
mid = aquarius::color_cnv[2],
high = aquarius::color_cnv[3],
midpoint = 0) +
ggplot2::theme_classic() +
ggplot2::theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
We represent a subset of genes on a heatmap :
features_oi = mark %>%
dplyr::filter(p_val_adj < 0.05) %>%
dplyr::filter(abs(avg_logFC) > 0.5) %>%
rownames()
length(features_oi)
## [1] 112
We prepare the scaled expression matrix :
mat_expression = Seurat::GetAssayData(subsobj, assay = "RNA", slot = "data")[features_oi, ]
mat_expression = Matrix::t(mat_expression)
mat_expression = dynutils::scale_quantile(mat_expression) # between 0 and 1
mat_expression = Matrix::t(mat_expression)
mat_expression = as.matrix(mat_expression) # not sparse
dim(mat_expression)
## [1] 112 2022
We prepare the heatmap annotation :
ha_top = ComplexHeatmap::HeatmapAnnotation(
cluster_type = subsobj$cluster_type,
sample_type = subsobj$sample_type,
cluster = subsobj$seurat_clusters,
col = list(cluster_type = color_markers,
sample_type = setNames(nm = c("HS", "HD"),
c("#C55F40", "#2C78E6")),
cluster = setNames(nm = levels(subsobj$seurat_clusters),
aquarius::gg_color_hue(length(levels(subsobj$seurat_clusters))))))
And the heatmap :
ht = ComplexHeatmap::Heatmap(mat_expression,
col = aquarius::color_cnv,
# Annotation
top_annotation = ha_top,
# Grouping
column_title = NULL,
cluster_rows = TRUE,
cluster_columns = TRUE,
show_column_names = FALSE,
# Visual aspect
show_heatmap_legend = TRUE,
border = TRUE)
ComplexHeatmap::draw(ht,
merge_legend = TRUE,
heatmap_legend_side = "bottom",
annotation_legend_side = "bottom")
We explore enrichment in gene sets for cluster 5.
genes_of_interest = rownames(mark_cluster5)
enrichr_results = aquarius::run_enrichr(gene_names = genes_of_interest,
gene_corresp = gene_corresp,
gene_sets = gene_sets[, c("gs_name", "ensembl_gene")],
make_plot = TRUE,
plot_title = "Up-regulated in cluster 5 compared to other ORS")
list_results[[group_name]]$enrichr_up = enrichr_results$ego
enrichr_results$plot +
ggplot2::theme(axis.text.y = element_text(size = 8))
We explore enrichment in gene sets for genes downregulated in cluster 5.
genes_of_interest = rownames(mark_not5)
enrichr_results = aquarius::run_enrichr(gene_names = genes_of_interest,
gene_corresp = gene_corresp,
gene_sets = gene_sets[, c("gs_name", "ensembl_gene")],
make_plot = TRUE,
plot_title = "Down-regulated in cluster 5 compared to other ORS")
list_results[[group_name]]$enrichr_dn = enrichr_results$ego
enrichr_results$plot +
ggplot2::theme(axis.text.y = element_text(size = 8))
We run a GSEA for all gene sets, from the full count matrix :
ranked_gene_list = aquarius::run_foldchange(Seurat::GetAssayData(subsobj, assay = "RNA", slot = "counts"),
group1 = colnames(subsobj)[subsobj@active.ident == 5],
group2 = colnames(subsobj)[subsobj@active.ident != 5])
names(ranked_gene_list) = gene_corresp$ID
gsea_results = aquarius::gsea_run(ranked_gene_list = ranked_gene_list,
gene_sets = gene_sets[, c("gs_name", "ensembl_gene")],
GSEA_p_val_thresh = 1)
list_results[[group_name]]$gsea = gsea_results
gsea_results@result %>%
dplyr::filter(pvalue < 0.05) %>%
dplyr::top_n(., n = 200, wt = abs(NES)) %>%
dplyr::mutate(too_long = ifelse(nchar(ID) > 60, yes = TRUE, no = FALSE)) %>%
dplyr::mutate(ID = stringr::str_sub(ID, end = 60)) %>%
dplyr::mutate(ID = ifelse(too_long, yes = paste0(ID, "..."), no = ID)) %>%
aquarius::gsea_plot(show_legend = TRUE) +
ggplot2::labs(title = "GSEA using all genes (count matrix)") +
ggplot2::theme(plot.title = element_text(size = 20))
We compute a fold change table to make a wordcloud :
fold_change = gene_corresp
colnames(fold_change) = c("gene_name", "ID")
rownames(fold_change) = fold_change$ID
fold_change$FC = ranked_gene_list[rownames(fold_change)]
fold_change$pct.1 = Seurat::GetAssayData(subsobj, assay = "RNA", slot = "counts")[, subsobj@active.ident == "5"] %>%
apply(., MARGIN = 1, FUN = function(one_row) {
return( mean(one_row != 0) )
})
fold_change$pct.2 = Seurat::GetAssayData(subsobj, assay = "RNA", slot = "counts")[, subsobj@active.ident != "5"] %>%
apply(., MARGIN = 1, FUN = function(one_row) {
return( mean(one_row != 0) )
})
fold_change$FC_x_pct = ifelse(fold_change$FC > 0,
yes = fold_change$FC * fold_change$pct.1,
no = fold_change$FC * fold_change$pct.2)
dim(fold_change) ; head(fold_change)
## [1] 16701 6
## gene_name ID FC pct.1 pct.2
## ENSG00000238009 AL627309.1 ENSG00000238009 0.998441346 0.009933775 0.004069767
## ENSG00000237491 AL669831.5 ENSG00000237491 0.985617306 0.162251656 0.055232558
## ENSG00000225880 LINC00115 ENSG00000225880 0.285723298 0.049668874 0.033720930
## ENSG00000230368 FAM41C ENSG00000230368 -1.335459391 0.036423841 0.065116279
## ENSG00000230699 AL645608.3 ENSG00000230699 -0.001558654 0.006622517 0.006395349
## ENSG00000187634 SAMD11 ENSG00000187634 -1.460990273 0.000000000 0.005232558
## FC_x_pct
## ENSG00000238009 9.918292e-03
## ENSG00000237491 1.599180e-01
## ENSG00000225880 1.419155e-02
## ENSG00000230368 -8.696015e-02
## ENSG00000230699 -9.968137e-06
## ENSG00000187634 -7.644717e-03
We make the gsea plot for some gene sets :
gs_oi = c("GOBP_CD4_POSITIVE_ALPHA_BETA_T_CELL_DIFFERENTIATION",
"GOBP_REGULATION_OF_ALPHA_BETA_T_CELL_DIFFERENTIATION",
"GOBP_ALPHA_BETA_T_CELL_DIFFERENTIATION",
"GOBP_POSITIVE_REGULATION_OF_ALPHA_BETA_T_CELL_ACTIVATION",
"GOBP_ALPHA_BETA_T_CELL_ACTIVATION")
plot_list = make_gsea_plot(gsea_results, gs_oi, fold_change, "FC_x_pct")
patchwork::wrap_plots(plot_list, ncol = 2, widths = c(2,1.5))
group_name = "IBL_HS_vs_HD"
We subset the Seurat object and change cell identities to sample type.
subsobj = subset(sobj, cluster_type == "IBL")
Seurat::Idents(subsobj) = subsobj$sample_type
table(subsobj$sample_type)
##
## HS HD
## 1284 226
We identify DE genes between HS and HD :
mark = Seurat::FindMarkers(subsobj, ident.1 = "HS", ident.2 = "HD")
mark = mark %>%
dplyr::filter(p_val_adj < 0.05) %>%
dplyr::arrange(-avg_logFC, pct.1 - pct.2)
list_results[[group_name]]$mark = mark
dim(mark)
## [1] 107 5
head(mark, n = 20)
## p_val avg_logFC pct.1 pct.2 p_val_adj
## LGALS7 1.216786e-49 1.6340915 0.685 0.190 2.032154e-45
## S100A7 2.913857e-11 1.5790394 0.289 0.080 4.866433e-07
## MIF 1.195470e-44 1.2201587 0.706 0.261 1.996554e-40
## LGALS7B 1.099708e-09 1.1642854 0.769 0.717 1.836623e-05
## MTRNR2L12 9.867130e-12 1.0580493 0.393 0.181 1.647909e-07
## MTRNR2L8 3.047443e-11 0.9124631 0.350 0.155 5.089535e-07
## MT-CO1 4.350400e-12 0.9105978 0.449 0.226 7.265603e-08
## FOSB 7.804218e-09 0.8760606 0.298 0.133 1.303382e-04
## NDUFA4L2 1.283179e-20 0.8201600 0.743 0.434 2.143037e-16
## MT-CO2 1.272279e-08 0.8121340 0.433 0.265 2.124832e-04
## S100A9 1.727393e-13 0.7487053 0.315 0.075 2.884919e-09
## KRT15 2.873193e-08 0.6857105 0.290 0.115 4.798519e-04
## CA2 1.116064e-16 0.6640743 0.460 0.173 1.863938e-12
## ARF5 3.584966e-27 0.6221800 0.433 0.044 5.987251e-23
## CKB 4.051432e-08 0.5928841 0.282 0.115 6.766297e-04
## MT-ND4 2.524958e-07 0.5633098 0.410 0.248 4.216932e-03
## MT-ATP6 1.537778e-06 0.5539256 0.496 0.358 2.568244e-02
## CST6 1.707459e-07 0.5234863 0.421 0.221 2.851627e-03
## IGKC 5.879378e-07 0.5045779 0.101 0.000 9.819150e-03
## MCL1 1.177137e-07 0.5038666 0.379 0.208 1.965937e-03
We explore enrichment in gene sets for HS population.
genes_of_interest = mark %>%
dplyr::filter(avg_logFC > 0) %>%
rownames()
enrichr_results = aquarius::run_enrichr(gene_names = genes_of_interest,
gene_corresp = gene_corresp,
gene_sets = gene_sets[, c("gs_name", "ensembl_gene")],
make_plot = TRUE,
plot_title = "Up-regulated in HS compared to HD")
list_results[[group_name]]$enrichr_hs = enrichr_results$ego
enrichr_results$plot +
ggplot2::theme(axis.text.y = element_text(size = 8))
We explore enrichment in gene sets for HD population.
genes_of_interest = mark %>%
dplyr::filter(avg_logFC < 0) %>%
rownames()
enrichr_results = aquarius::run_enrichr(gene_names = genes_of_interest,
gene_corresp = gene_corresp,
gene_sets = gene_sets[, c("gs_name", "ensembl_gene")],
make_plot = TRUE,
plot_title = "Down-regulated in HS compared to HD")
list_results[[group_name]]$enrichr_hd = enrichr_results$ego
enrichr_results$plot +
ggplot2::theme(axis.text.y = element_text(size = 8))
We run a GSEA for all gene sets, from the full count matrix :
ranked_gene_list = aquarius::run_foldchange(Seurat::GetAssayData(subsobj, assay = "RNA", slot = "counts"),
group1 = colnames(subsobj)[subsobj@active.ident %in% "HS"],
group2 = colnames(subsobj)[subsobj@active.ident %in% "HD"])
names(ranked_gene_list) = gene_corresp$ID
gsea_results = aquarius::gsea_run(ranked_gene_list = ranked_gene_list,
gene_sets = gene_sets[, c("gs_name", "ensembl_gene")],
GSEA_p_val_thresh = 1)
list_results[[group_name]]$gsea = gsea_results
gsea_results@result %>%
dplyr::filter(pvalue < 0.05) %>%
dplyr::top_n(., n = 200, wt = abs(NES)) %>%
dplyr::mutate(too_long = ifelse(nchar(ID) > 60, yes = TRUE, no = FALSE)) %>%
dplyr::mutate(ID = stringr::str_sub(ID, end = 60)) %>%
dplyr::mutate(ID = ifelse(too_long, yes = paste0(ID, "..."), no = ID)) %>%
aquarius::gsea_plot(show_legend = TRUE) +
ggplot2::labs(title = "GSEA using all genes (count matrix)") +
ggplot2::theme(plot.title = element_text(size = 20))
We compute a fold change table to make a wordcloud :
fold_change = gene_corresp
colnames(fold_change) = c("gene_name", "ID")
rownames(fold_change) = fold_change$ID
fold_change$FC = ranked_gene_list[rownames(fold_change)]
fold_change$pct.1 = Seurat::GetAssayData(subsobj, assay = "RNA", slot = "counts")[, subsobj@active.ident == "HS"] %>%
apply(., MARGIN = 1, FUN = function(one_row) {
return( mean(one_row != 0) )
})
fold_change$pct.2 = Seurat::GetAssayData(subsobj, assay = "RNA", slot = "counts")[, subsobj@active.ident == "HD"] %>%
apply(., MARGIN = 1, FUN = function(one_row) {
return( mean(one_row != 0) )
})
fold_change$FC_x_pct = ifelse(fold_change$FC > 0,
yes = fold_change$FC * fold_change$pct.1,
no = fold_change$FC * fold_change$pct.2)
dim(fold_change) ; head(fold_change)
## [1] 16701 6
## gene_name ID FC pct.1 pct.2
## ENSG00000238009 AL627309.1 ENSG00000238009 -0.79640583 0.0046728972 0.000000000
## ENSG00000237491 AL669831.5 ENSG00000237491 0.04968149 0.0919003115 0.044247788
## ENSG00000225880 LINC00115 ENSG00000225880 -1.01879825 0.0210280374 0.017699115
## ENSG00000230368 FAM41C ENSG00000230368 1.12415970 0.0389408100 0.004424779
## ENSG00000230699 AL645608.3 ENSG00000230699 -0.60376075 0.0054517134 0.000000000
## ENSG00000187634 SAMD11 ENSG00000187634 -3.60376075 0.0007788162 0.004424779
## FC_x_pct
## ENSG00000238009 0.000000000
## ENSG00000237491 0.004565744
## ENSG00000225880 -0.018031827
## ENSG00000230368 0.043775689
## ENSG00000230699 0.000000000
## ENSG00000187634 -0.015945844
We make the gsea plot for some gene sets :
gs_oi = c("REACTOME_EXPORT_OF_VIRAL_RIBONUCLEOPROTEINS_FROM_NUCLEUS",
"GOBP_MICROTUBULE_ANCHORING",
"GOBP_HISTONE_H3_K4_TRIMETHYLATION",
"GOBP_HISTONE_H3_K4_MONOMETHYLATION",
"GOBP_MITOPHAGY",
"REACTOME_SEROTONIN_NEUROTRANSMITTER_RELEASE_CYCLE",
"REACTOME_DOPAMINE_NEUROTRANSMITTER_RELEASE_CYCLE")
plot_list = make_gsea_plot(gsea_results, gs_oi, fold_change, "FC_x_pct")
patchwork::wrap_plots(plot_list, ncol = 2, widths = c(2,1.5))
group_name = "ORS_HS_vs_HD"
We subset the Seurat object and change cell identities to sample type.
subsobj = subset(sobj, cluster_type == "ORS")
Seurat::Idents(subsobj) = subsobj$sample_type
table(subsobj$sample_type)
##
## HS HD
## 1563 459
We identify DE genes between HS and HD :
mark = Seurat::FindMarkers(subsobj, ident.1 = "HS", ident.2 = "HD")
mark = mark %>%
dplyr::filter(p_val_adj < 0.05) %>%
dplyr::arrange(-avg_logFC, pct.1 - pct.2)
list_results[[group_name]]$mark = mark
dim(mark)
## [1] 145 5
head(mark, n = 20)
## p_val avg_logFC pct.1 pct.2 p_val_adj
## S100A7 1.177379e-34 1.3981213 0.390 0.083 1.966341e-30
## S100A9 8.042694e-80 1.3882495 0.831 0.346 1.343210e-75
## S100A8 6.069068e-63 1.2071418 0.784 0.353 1.013595e-58
## RPS26 8.289007e-157 1.0156704 0.971 0.950 1.384347e-152
## LGALS7B 1.751974e-26 0.9144582 0.636 0.503 2.925971e-22
## CCL2 4.001592e-35 0.8469394 0.592 0.279 6.683058e-31
## LTF 1.036558e-17 0.7369496 0.183 0.022 1.731156e-13
## MTRNR2L8 2.332458e-41 0.7044528 0.909 0.904 3.895438e-37
## IFI27 4.002962e-19 0.6856359 0.194 0.024 6.685347e-15
## TIMP1 1.883812e-25 0.6120579 0.566 0.342 3.146155e-21
## IFITM3 2.955422e-49 0.5618141 0.889 0.819 4.935850e-45
## CD74 9.221913e-27 0.5594045 0.347 0.096 1.540152e-22
## HLA-C 1.036948e-24 0.5483711 0.731 0.547 1.731807e-20
## MIF 1.151305e-28 0.5462406 0.403 0.159 1.922795e-24
## AKR1B10 7.040926e-44 0.5197609 0.372 0.033 1.175905e-39
## MTRNR2L12 1.018482e-25 0.4788943 0.951 0.943 1.700967e-21
## MTRNR2L10 1.055574e-11 0.4753293 0.589 0.538 1.762914e-07
## CXCL14 2.007187e-14 0.4495196 0.882 0.802 3.352203e-10
## LGALS7 5.171905e-12 0.4358938 0.232 0.096 8.637599e-08
## IL1R2 9.715791e-13 0.4197815 0.226 0.085 1.622634e-08
We explore enrichment in gene sets for HS population.
genes_of_interest = mark %>%
dplyr::filter(avg_logFC > 0) %>%
rownames()
enrichr_results = aquarius::run_enrichr(gene_names = genes_of_interest,
gene_corresp = gene_corresp,
gene_sets = gene_sets[, c("gs_name", "ensembl_gene")],
make_plot = TRUE,
plot_title = "Up-regulated in HS compared to HD")
list_results[[group_name]]$enrichr_hs = enrichr_results$ego
enrichr_results$plot +
ggplot2::theme(axis.text.y = element_text(size = 8))
We explore enrichment in gene sets for HD population.
genes_of_interest = mark %>%
dplyr::filter(avg_logFC < 0) %>%
rownames()
enrichr_results = aquarius::run_enrichr(gene_names = genes_of_interest,
gene_corresp = gene_corresp,
gene_sets = gene_sets[, c("gs_name", "ensembl_gene")],
make_plot = TRUE,
plot_title = "Down-regulated in HS compared to HD")
list_results[[group_name]]$enrichr_hd = enrichr_results$ego
enrichr_results$plot +
ggplot2::theme(axis.text.y = element_text(size = 8))
We run a GSEA for all gene sets, from the full count matrix :
ranked_gene_list = aquarius::run_foldchange(Seurat::GetAssayData(subsobj, assay = "RNA", slot = "counts"),
group1 = colnames(subsobj)[subsobj@active.ident %in% "HS"],
group2 = colnames(subsobj)[subsobj@active.ident %in% "HD"])
names(ranked_gene_list) = gene_corresp$ID
gsea_results = aquarius::gsea_run(ranked_gene_list = ranked_gene_list,
gene_sets = gene_sets[, c("gs_name", "ensembl_gene")],
GSEA_p_val_thresh = 1)
list_results[[group_name]]$gsea = gsea_results
gsea_results@result %>%
dplyr::filter(pvalue < 0.05) %>%
dplyr::top_n(., n = 200, wt = abs(NES)) %>%
dplyr::mutate(too_long = ifelse(nchar(ID) > 60, yes = TRUE, no = FALSE)) %>%
dplyr::mutate(ID = stringr::str_sub(ID, end = 60)) %>%
dplyr::mutate(ID = ifelse(too_long, yes = paste0(ID, "..."), no = ID)) %>%
aquarius::gsea_plot(show_legend = TRUE) +
ggplot2::labs(title = "GSEA using all genes (count matrix)") +
ggplot2::theme(plot.title = element_text(size = 20))
We compute a fold change table to make a wordcloud :
fold_change = gene_corresp
colnames(fold_change) = c("gene_name", "ID")
rownames(fold_change) = fold_change$ID
fold_change$FC = ranked_gene_list[rownames(fold_change)]
fold_change$pct.1 = Seurat::GetAssayData(subsobj, assay = "RNA", slot = "counts")[, subsobj@active.ident == "HS"] %>%
apply(., MARGIN = 1, FUN = function(one_row) {
return( mean(one_row != 0) )
})
fold_change$pct.2 = Seurat::GetAssayData(subsobj, assay = "RNA", slot = "counts")[, subsobj@active.ident == "HD"] %>%
apply(., MARGIN = 1, FUN = function(one_row) {
return( mean(one_row != 0) )
})
fold_change$FC_x_pct = ifelse(fold_change$FC > 0,
yes = fold_change$FC * fold_change$pct.1,
no = fold_change$FC * fold_change$pct.2)
dim(fold_change) ; head(fold_change)
## [1] 16701 6
## gene_name ID FC pct.1 pct.2
## ENSG00000238009 AL627309.1 ENSG00000238009 -0.2330146 0.005118362 0.004357298
## ENSG00000237491 AL669831.5 ENSG00000237491 0.7769695 0.079334613 0.043572985
## ENSG00000225880 LINC00115 ENSG00000225880 0.3519479 0.038387716 0.028322440
## ENSG00000230368 FAM41C ENSG00000230368 0.1550558 0.062699936 0.054466231
## ENSG00000230699 AL645608.3 ENSG00000230699 -0.8179771 0.005758157 0.008714597
## ENSG00000187634 SAMD11 ENSG00000187634 -1.3325503 0.003198976 0.008714597
## FC_x_pct
## ENSG00000238009 -0.001015314
## ENSG00000237491 0.061640571
## ENSG00000225880 0.013510475
## ENSG00000230368 0.009721990
## ENSG00000230699 -0.007128341
## ENSG00000187634 -0.011612639
We make the gsea plot for some gene sets :
gs_oi = c("GOBP_RESPONSE_TO_MOLECULE_OF_BACTERIAL_ORIGIN",
"GOBP_HUMORAL_IMMUNE_RESPONSE",
"HALLMARK_APOPTOSIS",
"GOBP_GRANULOCYTE_CHEMOTAXIS",
"GOBP_GRANULOCYTE_MIGRATION",
"HALLMARK_INTERFERON_ALPHA_RESPONSE",
"HALLMARK_INTERFERON_GAMMA_RESPONSE",
"HALLMARK_TNFA_SIGNALING_VIA_NFKB",
"HALLMARK_HYPOXIA",
"GOBP_RESPONSE_TO_TRANSFORMING_GROWTH_FACTOR_BETA")
plot_list = make_gsea_plot(gsea_results, gs_oi, fold_change, "FC_x_pct")
patchwork::wrap_plots(plot_list, ncol = 2, widths = c(2,1.5))
We represent all differentially expressed genes on a heatmap. First, we extract all DE genes :
features_oi = c(mark_to_label$gene_name,
rownames(list_results$IBL_HS_vs_HD$mark),
rownames(list_results$ORS_HS_vs_HD$mark)) %>%
unique()
length(features_oi)
## [1] 259
We prepare the scaled expression matrix :
mat_expression = Seurat::GetAssayData(sobj, assay = "RNA", slot = "data")[features_oi, ]
mat_expression = Matrix::t(mat_expression)
mat_expression = dynutils::scale_quantile(mat_expression) # between 0 and 1
mat_expression = Matrix::t(mat_expression)
mat_expression = as.matrix(mat_expression) # not sparse
dim(mat_expression)
## [1] 259 3532
We prepare the heatmap annotation :
ha_top = ComplexHeatmap::HeatmapAnnotation(
cluster_type = sobj$cluster_type,
sample_type = sobj$sample_type,
cluster = sobj$seurat_clusters,
col = list(cluster_type = color_markers,
sample_type = setNames(nm = c("HS", "HD"),
c("#C55F40", "#2C78E6")),
cluster = setNames(nm = levels(sobj$seurat_clusters),
aquarius::gg_color_hue(length(levels(sobj$seurat_clusters))))))
And the heatmap :
sobj$cell_group = paste0(sobj$cluster_type, sobj$sample_type) %>%
as.factor()
ht = ComplexHeatmap::Heatmap(mat_expression,
col = aquarius::color_cnv,
# Annotation
top_annotation = ha_top,
# Grouping
column_order = sobj@meta.data %>%
dplyr::arrange(cluster_type, sample_type, seurat_clusters) %>%
rownames(),
column_split = sobj$cell_group,
column_gap = unit(c(0.01, 2, 0.01), "mm"),
column_title = NULL,
cluster_rows = TRUE,
cluster_columns = FALSE,
show_column_names = FALSE,
# Visual aspect
show_heatmap_legend = TRUE,
border = TRUE)
ComplexHeatmap::draw(ht,
merge_legend = TRUE,
heatmap_legend_side = "bottom",
annotation_legend_side = "bottom")
We save the list of results :
saveRDS(list_results, file = paste0(out_dir, "/", save_name, "_list_results.rds"))
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.6 LTS
##
## Matrix products: default
## BLAS: /usr/local/lib/R/lib/libRblas.so
## LAPACK: /usr/local/lib/R/lib/libRlapack.so
##
## locale:
## [1] C
##
## attached base packages:
## [1] parallel stats4 grid stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] org.Mm.eg.db_3.10.0 AnnotationDbi_1.48.0 IRanges_2.20.2
## [4] S4Vectors_0.24.4 Biobase_2.46.0 BiocGenerics_0.32.0
## [7] ComplexHeatmap_2.14.0 ggplot2_3.3.5 patchwork_1.1.2
## [10] dplyr_1.0.7
##
## loaded via a namespace (and not attached):
## [1] softImpute_1.4 graphlayouts_0.7.0
## [3] pbapply_1.4-2 lattice_0.20-41
## [5] haven_2.3.1 vctrs_0.3.8
## [7] usethis_2.0.1 dynwrap_1.2.1
## [9] blob_1.2.1 survival_3.2-13
## [11] prodlim_2019.11.13 dynutils_1.0.5
## [13] later_1.3.0 DBI_1.1.1
## [15] R.utils_2.11.0 SingleCellExperiment_1.8.0
## [17] rappdirs_0.3.3 uwot_0.1.8
## [19] dqrng_0.2.1 jpeg_0.1-8.1
## [21] zlibbioc_1.32.0 pspline_1.0-18
## [23] pcaMethods_1.78.0 mvtnorm_1.1-1
## [25] htmlwidgets_1.5.4 GlobalOptions_0.1.2
## [27] future_1.22.1 UpSetR_1.4.0
## [29] laeken_0.5.2 leiden_0.3.3
## [31] clustree_0.4.3 scater_1.14.6
## [33] irlba_2.3.3 markdown_1.1
## [35] DEoptimR_1.0-9 tidygraph_1.1.2
## [37] Rcpp_1.0.9 readr_2.0.2
## [39] KernSmooth_2.23-17 carrier_0.1.0
## [41] promises_1.1.0 gdata_2.18.0
## [43] DelayedArray_0.12.3 limma_3.42.2
## [45] graph_1.64.0 RcppParallel_5.1.4
## [47] Hmisc_4.4-0 fs_1.5.2
## [49] RSpectra_0.16-0 fastmatch_1.1-0
## [51] ranger_0.12.1 digest_0.6.25
## [53] png_0.1-7 sctransform_0.2.1
## [55] cowplot_1.0.0 DOSE_3.12.0
## [57] here_1.0.1 TInGa_0.0.0.9000
## [59] ggraph_2.0.3 pkgconfig_2.0.3
## [61] GO.db_3.10.0 DelayedMatrixStats_1.8.0
## [63] gower_0.2.1 ggbeeswarm_0.6.0
## [65] iterators_1.0.12 DropletUtils_1.6.1
## [67] reticulate_1.26 clusterProfiler_3.14.3
## [69] SummarizedExperiment_1.16.1 circlize_0.4.15
## [71] beeswarm_0.4.0 GetoptLong_1.0.5
## [73] xfun_0.35 bslib_0.3.1
## [75] zoo_1.8-10 tidyselect_1.1.0
## [77] reshape2_1.4.4 purrr_0.3.4
## [79] ica_1.0-2 pcaPP_1.9-73
## [81] viridisLite_0.3.0 rtracklayer_1.46.0
## [83] rlang_1.0.2 hexbin_1.28.1
## [85] jquerylib_0.1.4 dyneval_0.9.9
## [87] glue_1.4.2 RColorBrewer_1.1-2
## [89] matrixStats_0.56.0 stringr_1.4.0
## [91] lava_1.6.7 europepmc_0.3
## [93] DESeq2_1.26.0 recipes_0.1.17
## [95] labeling_0.3 httpuv_1.5.2
## [97] class_7.3-17 BiocNeighbors_1.4.2
## [99] DO.db_2.9 annotate_1.64.0
## [101] jsonlite_1.7.2 XVector_0.26.0
## [103] bit_4.0.4 mime_0.9
## [105] aquarius_0.1.5 Rsamtools_2.2.3
## [107] gridExtra_2.3 gplots_3.0.3
## [109] stringi_1.4.6 processx_3.5.2
## [111] gsl_2.1-6 bitops_1.0-6
## [113] cli_3.0.1 batchelor_1.2.4
## [115] RSQLite_2.2.0 randomForest_4.6-14
## [117] tidyr_1.1.4 data.table_1.14.2
## [119] rstudioapi_0.13 GenomicAlignments_1.22.1
## [121] nlme_3.1-147 qvalue_2.18.0
## [123] scran_1.14.6 locfit_1.5-9.4
## [125] scDblFinder_1.1.8 listenv_0.8.0
## [127] ggthemes_4.2.4 gridGraphics_0.5-0
## [129] R.oo_1.24.0 dbplyr_1.4.4
## [131] TTR_0.24.2 readxl_1.3.1
## [133] lifecycle_1.0.1 timeDate_3043.102
## [135] ggpattern_0.3.1 munsell_0.5.0
## [137] cellranger_1.1.0 R.methodsS3_1.8.1
## [139] proxyC_0.1.5 visNetwork_2.0.9
## [141] caTools_1.18.0 codetools_0.2-16
## [143] ggwordcloud_0.5.0 GenomeInfoDb_1.22.1
## [145] vipor_0.4.5 lmtest_0.9-38
## [147] msigdbr_7.5.1 htmlTable_1.13.3
## [149] triebeard_0.3.0 lsei_1.2-0
## [151] xtable_1.8-4 ROCR_1.0-7
## [153] BiocManager_1.30.10 scatterplot3d_0.3-41
## [155] abind_1.4-5 farver_2.0.3
## [157] parallelly_1.28.1 RANN_2.6.1
## [159] askpass_1.1 GenomicRanges_1.38.0
## [161] RcppAnnoy_0.0.16 tibble_3.1.5
## [163] ggdendro_0.1-20 cluster_2.1.0
## [165] future.apply_1.5.0 Seurat_3.1.5
## [167] dendextend_1.15.1 Matrix_1.3-2
## [169] ellipsis_0.3.2 prettyunits_1.1.1
## [171] lubridate_1.7.9 ggridges_0.5.2
## [173] igraph_1.2.5 RcppEigen_0.3.3.7.0
## [175] fgsea_1.12.0 remotes_2.4.2
## [177] scBFA_1.0.0 destiny_3.0.1
## [179] VIM_6.1.1 testthat_3.1.0
## [181] htmltools_0.5.2 BiocFileCache_1.10.2
## [183] yaml_2.2.1 utf8_1.1.4
## [185] plotly_4.9.2.1 XML_3.99-0.3
## [187] ModelMetrics_1.2.2.2 e1071_1.7-3
## [189] foreign_0.8-76 withr_2.5.0
## [191] fitdistrplus_1.0-14 BiocParallel_1.20.1
## [193] xgboost_1.4.1.1 bit64_4.0.5
## [195] foreach_1.5.0 robustbase_0.93-9
## [197] Biostrings_2.54.0 GOSemSim_2.13.1
## [199] rsvd_1.0.3 memoise_2.0.0
## [201] evaluate_0.18 forcats_0.5.0
## [203] rio_0.5.16 geneplotter_1.64.0
## [205] tzdb_0.1.2 caret_6.0-86
## [207] ps_1.6.0 DiagrammeR_1.0.6.1
## [209] curl_4.3 fdrtool_1.2.15
## [211] fansi_0.4.1 highr_0.8
## [213] urltools_1.7.3 xts_0.12.1
## [215] GSEABase_1.48.0 acepack_1.4.1
## [217] edgeR_3.28.1 checkmate_2.0.0
## [219] scds_1.2.0 cachem_1.0.6
## [221] npsurv_0.4-0 babelgene_22.3
## [223] rjson_0.2.20 openxlsx_4.1.5
## [225] ggrepel_0.9.1 clue_0.3-60
## [227] rprojroot_2.0.2 stabledist_0.7-1
## [229] tools_3.6.3 sass_0.4.0
## [231] nichenetr_1.1.1 magrittr_2.0.1
## [233] RCurl_1.98-1.2 proxy_0.4-24
## [235] car_3.0-11 ape_5.3
## [237] ggplotify_0.0.5 xml2_1.3.2
## [239] httr_1.4.2 assertthat_0.2.1
## [241] rmarkdown_2.18 boot_1.3-25
## [243] globals_0.14.0 R6_2.4.1
## [245] Rhdf5lib_1.8.0 nnet_7.3-14
## [247] RcppHNSW_0.2.0 progress_1.2.2
## [249] genefilter_1.68.0 statmod_1.4.34
## [251] gtools_3.8.2 shape_1.4.6
## [253] HDF5Array_1.14.4 BiocSingular_1.2.2
## [255] rhdf5_2.30.1 splines_3.6.3
## [257] AUCell_1.8.0 carData_3.0-4
## [259] colorspace_1.4-1 generics_0.1.0
## [261] base64enc_0.1-3 dynfeature_1.0.0
## [263] smoother_1.1 gridtext_0.1.1
## [265] pillar_1.6.3 tweenr_1.0.1
## [267] sp_1.4-1 ggplot.multistats_1.0.0
## [269] rvcheck_0.1.8 GenomeInfoDbData_1.2.2
## [271] plyr_1.8.6 gtable_0.3.0
## [273] zip_2.2.0 knitr_1.41
## [275] latticeExtra_0.6-29 biomaRt_2.42.1
## [277] fastmap_1.1.0 ADGofTest_0.3
## [279] copula_1.0-0 doParallel_1.0.15
## [281] vcd_1.4-8 babelwhale_1.0.1
## [283] openssl_1.4.1 scales_1.1.1
## [285] backports_1.2.1 ipred_0.9-12
## [287] enrichplot_1.6.1 hms_1.1.1
## [289] ggforce_0.3.1 Rtsne_0.15
## [291] shiny_1.7.1 numDeriv_2016.8-1.1
## [293] polyclip_1.10-0 lazyeval_0.2.2
## [295] Formula_1.2-3 tsne_0.1-3
## [297] crayon_1.3.4 MASS_7.3-54
## [299] pROC_1.16.2 viridis_0.5.1
## [301] dynparam_1.0.0 rpart_4.1-15
## [303] zinbwave_1.8.0 compiler_3.6.3
## [305] ggtext_0.1.0